Divergent selection predating the Last Glacial Maximum mainly acted on macro‐phenotypes in Norway spruce

Abstract The current distribution and population structure of many species were, to a large extent, shaped by cycles of isolation in glacial refugia and subsequent population expansions. Isolation in and postglacial expansion through heterogeneous environments led to either neutral or adaptive divergence. Norway spruce is no exception, and its current distribution is the consequence of a constant interplay between evolutionary and demographic processes. We investigated population differentiation and adaptation of Norway spruce for juvenile growth, diameter of the stem, wood density, and tracheid traits at breast height. Data from 4461 phenotyped and genotyped Norway spruce from 396 half‐sib families in two progeny tests were used to test for divergent selection in the framework of Q ST vs. F ST. We show that the macroscopic resultant trait (stem diameter), unlike its microscopic components (tracheid dimensions) and juvenile growth, was under divergent selection that predated the Last Glacial Maximum. Altogether, the current variation in these phenotypic traits in Norway spruce is better explained by local adaptation to ancestral environments than to current ones, where populations were partly preadapted, mainly through growth‐related traits.

and selection is known to play an important role in shaping species distribution (Case et al., 2005;Li et al., 2022;Savolainen et al., 2007;Wisz et al., 2013). Among others, wood is one of the main tissues undergoing selection since radial growth of the woody stem provides the mechanical support for woody plants to grow tall and gain a competitive advantage (Falster & Westoby, 2003Givnish, 1982;Iwasa et al., 1985). However, adaptation for upright growth is often associated with modified wood properties and evolutionary tradeoffs (Chave et al., 2009;Poorter et al., 2006). Wood density is a good example of such a growth-survival trade-off (Kitajima, 1994;Poorter et al., 2010), where low wood density, resulting from wider lumina and, in hardwood species, wider vessels (e.g., Coomes et al., 2008;Enquist et al., 1999;Preston et al., 2006), is associated with fast growth. Conversely, high wood density provides higher biomechanical safety (e.g., Niklas, 1993;Van Gelder et al., 2006) and resistance against pathogens (Augspurger & Kelly, 1984), and is thus associated with higher survival rates. The growth-survival trade-off is a common ecological trade-off observed in both tropical (e.g., Falster & Westoby, 2005;Poorter, 2008) and boreal species (e.g., Pretzsch et al., 2018;Rossi et al., 2015).
Norway spruce (Picea abies (L.) H. Karst.) is a dominant boreal conifer species that thrives nowadays in a vast northern domain spanning from Norway to central Russia, as well as in smaller European domains, such as the Alps and the Carpathians (Li et al., 2022;Nota et al., 2022;Tsuda et al., 2016;Vendramin et al., 2000). Repeated demographic movements associated with climatic cycles and ensuing secondary contacts led to the division of Norway spruce into seven genetic clusters : Northern Fennoscandia (NFE), Central and South Sweden (CSE), Russia-Baltics (RBA), Northern Poland (NPL), Central Europe (CEU), Alpine (ALP), and Carpathian (ROM) domains. In Southern Sweden, recent hybrids between CSE and ALP clusters can be encountered (CSE-ALP). The timing of the divergence was given by that of the population structure that long predated the Last Glacial Maximum (10-15 million years ago ;Chen et al., 2019;Li et al., 2022). The demographic history leading to this population structure went along a history of successive adaptations to biotic and abiotic factors. Following the hypothesis that most adaptations are reflected by cluster-specific allele frequencies , the population structure of Norway spruce should partly be associated with signatures of past adaptations.
Investigating this association is essential to understand how much of the current phenotypic variation is due to adaptations predating the Last Glacial Maximum.
The comparison of neutral (F ST ) and quantitative trait (Q ST ) genetic divergence offers a theoretically well-grounded, accessible, and powerful tool for assessing the role of natural selection in such a population differentiation scenario (Lamy et al., 2012).
Under neutrality, it is expected that F ST = Q ST (Whitlock, 2008), even if this framework needs to be updated to overcome some limitations (e.g., De Villemereuil et al., 2020). By implementing a robust way of estimating F ST and Q ST , we used this framework to assess whether a phenotype diverges more between the genetic clusters than expected under random genetic drift. Additionally, in order to account for possible evolutionary trade-offs between wood-related traits, we analyzed several phenotypes: stem diameter and ring widths, wood density, and tracheid traits defining it.
Since the nature of selection between trees is expected to change over ontogenetic stages (Miriti, 2006), we studied the development of these phenotypes with cambial age/ring number from increment cores sampled from two progeny trials of the Swedish breeding program. We have conducted the analyses on a large dataset of 4461 20-year-old trees from 396 half-sib families, and 118,145 SNPs from exome sequencing. We show that the current phenotypic variation is partly explained by divergent selection that predated the Last Glacial Maximum, which mainly acted on macro-phenotypes (stem diameter) rather than micro-phenotypes (tracheid dimensions).

| Plant material
The study is based on plant material first presented by Chen et al. (2014), which is here analyzed with special attention to population structure and selection. The material originates from two prog-

| Population structure
Population structure was previously assessed at a fine scale on 1672 individuals using 399,801 unlinked putatively noncoding SNPs (see details in Chen et al., 2019), categorizing families into the following genetic clusters (ordered by decreasing latitude): Central and South Sweden (CSE), Russia-Baltics (RBA), Northern Poland (NPL), Central Europe (CEU), Alps (ALP), Carpathians (ROM), and hybrids between CSE and ALP (CSE-ALP) (for tree geographical origins, see figure 1 in Milesi et al., 2019). The families examined in this study followed the same genetic clustering (Figure 1a,b). Individuals from the ROM genetic cluster were removed, as our sample included only three such families, whereas all other clusters included at least 10 families (Table S1).

| Phenotype measurements
We reused the phenotypic dataset presented by Chen et al. (2014).
Briefly, increment cores of 12 mm diameter, extracted from the northern side of the stem, were collected in 2010 and 2011 at breast height (BH, 1.3 m above ground) from six trees per halfsib family per trial. The usual cutting age is 60-70 years or more, while the cores were extracted from 20-or 21-year-old trees: only juvenile stages of wood characteristics were then captured by the cores. The increment cores were analyzed by SilviScan (Evans, 1994(Evans, , 2006 at Innventia, now RISE (Stockholm, Sweden), assessing high-resolution pith-to-bark radial variation in growth, tracheid, and wood properties. In this study, we focused on the development with cambial age (ring number) of ring width (RW, in mm), and averages per ring of wood density (WD, in kg/m 3 ), tracheid radial width (TRW, in μm), tracheid tangential width (TTW, in μm), tracheid wall thickness (TWT, in μm), diameter at breast height (DBH, in mm) calculated from ring widths, and the age when the trees reached breast height (ABH). Although ABH is related to upright growth, it can be estimated accurately with high-quality increment cores from old trees. The measurements, identification of annual rings, and the importance of ABH are discussed in Lundqvist et al. (2018).
The two innermost rings were removed as the extremely curved structure of the xylem does not allow the same high data quality as for rings further out. In addition, only data up to the 13th ring was kept for analyses since only fast-growing offspring with low ABH are represented at the highest cambial ages. To avoid the effect of the thinning performed, we removed all annual rings formed after 2008.
The study then covered six traits (RW, WD, TRW, TTW, TWT, and DBH) for 11 rings per offspring in most of the cases and their ABH.

| Exome capture and SNP filtering
We reused the genotypic information of the plant material available on BioProject PRJNA511374  and BioProject PRJNA731384 (Chen et al., 2021). The genotypic dataset is the exome capture (target resequencing) of the mother of each half-sib family, where 120-bp-length probes are designed to capture 26,219 genes of Picea abies (Vidalis et al., 2018). DNA extraction protocol of the 518 maternal trees is detailed in Chen et al. (2021), and SNP calling is detailed in Li et al. (2022).
As a complementary filter and as an alternative way to account for minor allele frequency (MAF), we additionally filtered out SNPs with less than five individuals per genotype to avoid a leverage effect. We performed a loose filtering for variant call rate (or missingness) with PLINK v1.9 (Chang et al., 2015), by removing individuals with more than 70% SNP missingness, resulting in 472 half-sib families, and by removing SNPs with more than 70% individual missingness, resulting in 118,145 SNPs (each mother is then genotyped

| Statistical model
The traits (DBH, WD, RW, TRW, TTW, and TWT) were normalized by their respective overall means (across all trees and rings 3-13) to allow comparison between traits. Then, they were analyzed at each cambial age separately, resulting in 6 traits over 11 rings, a total of 66 traits in addition to ABH.
To study the family effect underlying the traits, we used the following mixed model: where "Y ijkl " is the trait of the l-th offspring of the k-th family in the j-th block at the i-th site, "μ" is the grand mean, "s i " is the fixed effect of the i-th site, "B ij " and "F k " are uncorrelated random effects of the j-th block nested in i-th site (with variance 2

B
) and the k-th half-sib family (with variance 2 F ), respectively, and "ε ijkl " is the residual fitted to a Gaussian distribution of variance σ 2 . Heritability was estimated as an intraclass correlation estimation of half-sib families (p.166, Falconer, 1989).
To investigate the effect of genetic clusters on traits jointly with the effect of ABH, we used the following linear mixed model: where "c l " is the fixed effect of the l-th offspring's ABH, "(c × g) l " is the fixed interaction effect between ABH and the genetic cluster of the l-th offspring, "(c × s) l : is the fixed interaction effect between ABH and the site of the l-th offspring, and the other terms are as defined above in model (1). Model (2) was applied to the 66 traits, as previously mentioned. Since the interaction between ABH and the site effect was always significant, we split the dataset per site, as recommended by Crawley (2007, pp. 323-386). The results are focusing on the Höreda trial, but the conclusions were the same in the Erikstorp trial.
Finally, to disentangle the effect of divergent selection and neutral divergence, we estimated Q ST as defined by Spitze (1993) with the following mixed model: Where "G k " is the random effect of the k-th family's genetic cluster All computations were done using R v4.0.4 (R Core Team, 2021), and all models were fitted with the R package lme4 (function lmer; Bates et al., 2015). Statistical significance of the fixed effects was assessed using a type II Wald chi-square test (function Anova of the R package car), and significance of a slope was assessed using a Student's t-test (function t.test of the R package stats). To account for the unbalance caused by consecutive filtering, we only kept families with at least four offspring per site. Applying all the previous filters, we ended up with 396 half-sib families and 4461 offspring per trait.
Multiple testing across the 66 traits were accounted for by considering a Bonferroni adjusted cut-off of p = 7.6 × 10 −4 .

| The genetic cluster effect was significant at all ages
At all cambial ages, phenotypic variation of all six phenotypes (DBH, WD, RW, TRW, TTW, and TWT) was associated with the geographical origin of the genetic cluster, as the genetic cluster effect was significant for all rings (χ 2 > 11.6, degree of freedom or df = 5, p = 0.041), except for the third ring TWT (χ 2 = 10.9, df = 5, p = 0.053). Every trait of the southern populations was consistent with a higher investment in growth than in wood density (higher DBH, RW, TRW, and TTW; lower WD and TWT).
Phenotypic variation exhibited a strong temporal pattern, an ontogenetic shift (Figure 1c). For all rings, phenotypic variation also depended on the geographical origin of the trees. The ontogeny was quantitatively measurable, as it followed the second axis of the PCA on the six phenotypes (41.2% by PC1 and 38.3% by PC2). As expected, the ontogeny was mainly driven by that of DBH, which alone explained 34% of PC2, and on average increased fourfold over the cambial ages investigated, while other phenotypes varied at most by 54% (for further details, see Lundqvist et al., 2018). As cambial age increased, ontogenetic shift was instead driven by other phenotypes, as the coefficient of variation of DBH was the only one decreasing over time (Student's t = −6.13, p < 0.001; Figure S1).
Overall, the ontogenetic shift resulted from a complex interplay of multiple phenotypes, all with heritabilities higher than 30% (Table 1), in accordance with previously reported values (Chen et al., 2014(Chen et al., , 2021Hannrup et al., 2004).

| Accounting for early growth rate ruled out environmental causes from the analysis of ontogeny
By definition, ontogenic tempo is given by the ABH. It had a relatively higher environmental variance, and therefore a relatively lower heritability of 17.31% than other traits (Table 1). The lower heritability was not due to a lower genetic variance, in fact, quite the opposite. ABH also depended on the geographical origin of the genetic cluster as the genetic cluster effect was significant (χ 2 = 102.8, df = 1, p < 10 −5 in model (1)). As mentioned earlier, ABH was latitudinally structured, half-sib families of southern origin having low ABH reflecting higher early growth rates.
Despite evidence supporting the genetic control of ABH, the genetic signals in the model (2) seem to have been absorbed by the genetic cluster and the family effects, and let the feature "ABH" only reflect environmental factors out of a lesser correlation with the response variable: when correcting for multiple testing, in model (2), the effects of ABH and the genetic clusters were significant for more than 70% of the phenotypes (47 of 66 for both), whereas their interaction was significant for <2% of the phenotypes (1 of 66; Table S2). In other words, ABH explained a large part of the phenotypic variation, while being orthogonal to the family or genetic cluster effect -the variable "ABH" was then likely only capturing the environmental part of ABH. In this particular statistical context, we will assume that ABH reflects environmental factors -else than the block effect already accounted for -impacting early growth and orthogonal to genetic factors (e.g., micro-environmental heterogeneity of the field, supposedly orthogonal to genetic factors because of the random block design). Under this assumption, it is possible to disentangle the effect of genetic factors on the phenotype of a "mature" tree from environmental factors giving advantages in early stages of life (i.e., acclimation). We will come back to this interpretation in the Discussion. On later rings, trees with a low ABH, that is, fast growers, ended up with a markedly larger DBH, lower WD and TWT, and slightly higher TRW and TTW. Fast growers had increasingly higher DBH over time, as the effect of ABH decreased from the 3rd ring to the 13th ring (from 5.36 × 10 −3 to −2.87 × 10 −2 ). On the opposite, fast growers had increasingly lower WD, as the effect of ABH increased from the 3rd ring to the 13th ring (from 4.93 × 10 −3 to 3.1 × 10 −2 ). Fast growers also had significantly larger 9th rings, as the effect of ABH on RW followed a convex curve matching its concave significance, reaching its largest negative value at the 9th ring (−5.82 × 10 −1 ). At all cambial ages, the effect of ABH was smaller for tracheid properties (TRW, TTW, and TWT) than for macro-phenotypes (DBH, WD) ( Figure S2). These results show that environmental factors favoring early growth had an impact on all later-stage phenotypes (except for RW), albeit stronger on macro-phenotypes (DBH, WD) than on micro-phenotypes (TWT, TRW, and TTW).

| Macro-phenotypes were also the ones showing signature of past divergent selection
In addition to all phenotypes being, as shown above, significantly as- (17.0% of variance explained) was mostly driven by TTW and TWT.
That the traits are structured by both family and genetic cluster can also be seen in the fact that intra-family (or intra-population) variance was significantly increased when the family (or population) structure was shuffled (see Supplement method C). Not only did phenotypic variation reflect the underlying population structure but it was also quantitatively associated with latitude.
The quantitative association between environmental gradients and phenotypic variation was partly due to selection (Figure 2c)

| DISCUSS ION
The Swedish breeding program is a mixture of local and recently introduced trees from other European countries. As such, although incomplete, the progeny trials reflect the population structure due to glaciation events and constitute an important source of information for evolutionary biology Li et al., 2022;Milesi et al., 2019). The breeding conditions constituted a unique opportunity to accurately correct phenotypes for environmental factors, by ruling out recruitment selection and considerably lowering interspecific competition. Furthermore, by correcting for population structure and using common garden experiments, we can rule out neutral and plastic processes so that the observed phenotypic differences can be interpreted as being the outcome of adaptation de Villemereuil et al. et al., 2020).
In our common garden experiments, growth exhibited a high level of population differentiation, as expected Savolainen et al., 2007). Population differentiation followed a latitudinal pattern, as previously reported Montwé et al., 2016;Rossi et al., 2015). The growth rate was generally differentiated among individual trees at early stages, supporting the fact that young ages are the most dynamic phases (Lundqvist et al., 2018).
Differences in apparent growth rate were mainly due to a longer period of growth rather than a faster daily growth (see discussion in Milesi et al., 2019), and despite the low heritability of growth rate, its high genetic and environmental variances suggest that it could be related to fitness (Houle, 1992;Caballero, 2020, p.235). This is also supported by the fact that the effect of ABH (here, only capturing the micro-environmental effect) was strongest on phenotypes on which signature of selection was also the strongest. Early advantages, possibly acquired by chance, seem to be re-invested in macrophenotypes, giving thus competitive advantages.

| The quantitative association between phenotypes and population structure predates the last glacial maximum
Both genetic and phenotypic variation had a clear latitudinal component, as previously reported in Norway spruce (e.g., Milesi et al., 2019) and in other species (e.g., Csilléry et al., 2020;Sang et al., 2019). Since local adaptation matches the environment  and postglacial recolonization routes followed environmental gradients (Gaggiotti et al., 2009;Li et al., 2022), population structure of Norway spruce is, as a result, strongly confounded with local adaptation. Population structure is thus a good predictor of selection in Norway spruce, as previously reported (Collignon et al., 2002;Kremer et al., 2010;Savolainen et al., 2007). However, local adaptation is not necessarily limited to recent population history. Some of the adaptations displayed by Norway spruce were associated with ancestral environments and some with current environments. If we consider that adaptation can be dated by the population structure it is associated with, knowing that the population structure of Norway spruce largely predates the Last Glacial Maximum (10-15 million years ago; Chen et al., 2019;Li et al., 2022), the observed phenotypic divergence is partly due to adaptations also predating the Last Glacial Maximum. We posit that the quantitative association among phenotype, latitude, and population structure would not be observed if there were no local adaptation to ancestral environments.
The last postglacial expansion of Norway spruce was thus mainly shaped by adaptations formerly acquired -not newly acquired. In other words, the phenotypic variation within the refugia explains in part the latitudinal variation observed nowadays at large geographical scale.
Because of the association between phenotypic variation and demographic processes, assessing the genetic control underlying adaptation of Norway spruce is challenging: not correcting for population structure in GWASs will lead to a high false-positive rate (Barton et al., 2019;Lawson et al., 2019), and correcting for population structure may remove true association, that is, lead to a high false-negative rate. A comprehensive assessment of the genetic variants underlying phenotypic variation is then not possible when there is such a strong latitudinal component. Increasing the number of samples or the number of SNPs remains inconclusive (Chen et al., 2021). A possible workaround suggested by Milesi et al. (2019) is to consider phenotypes that are less latitudinally associated, such as the longitudinally associated bud burst, and correct other phenotypes by the association with population structure thereby inferred.
Alternatively, we can consider several species jointly as suggested by Alberto et al. (2013).

| Divergent selection mostly affected resultant traits, not components
The strongest genetic divergence was observed for growth-related traits, and secondly for wood density, suggesting that most of the

| CON CLUS ION
As expected, phenotypic variation reflected the population structure of Norway spruce. The fact that the phenotypes were latitudinally ordered suggests that the origin of population was strongly confounded with adaptation. Teasing apart the demographic history and the adaptation history is challenging, more so that it has been demonstrated that population history had a strong impact on phenotypic divergence in Norway spruce (Lagercrantz & Ryman, 1990;Milesi et al., 2019). Divergence of phenotypic traits might in part be the results of past adaptations, which can be dated by the population structure it is associated with. In accordance with previous studies (Martínez-Sancho et al., 2021;Milesi et al., 2019;Savolainen et al., 2013), our results support the hypothesis that the current distribution of Norway spruce reflects in part adaptations to ancestral environments, as some signatures of selection were associated with population structure that predated the Last Glacial Maximum. The divergent selection affected more growth-related macro-phenotypes than micro-phenotypes at tracheid level -presumably invisible to selection.

ACK N OWLED G M ENTS
The authors acknowledge the Swedish research program Bio4Energy and Skogforsk for access to phenotypic data, and the Swedish

CO N FLI C T O F I NTE R E S T
The authors declare that there is no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The exome capture raw reads and the RNA-seq data have been deposited in NCBI's sequence read archive (SRA) under accession number PRJNA511374  and PRJNA731384 (Chen et al., 2021). For question on availability of the raw phenotypic data, contact RGG. Maria Rosario García-Gil https://orcid.org/0000-0002-6834-6708